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ABSTRACT 


An  algorithm  for  multiple  target  tracking  and  data 
correlation  is  described.  A  general  description  of  the 
problem  and  solution  is  first  given.  More  specific 
discusssions  on  tracking  with  a  passive  infrared  sensor  then 
follow.  An  example  is  presented  to  illustrate  the  trade-off 
between  algorithm  complexity,  performance,  and  processing 
requirements . 
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INTRODUCTION 


I . 

Multiple  target  tracking  is  a  classical  problem  with 
both  civilian  and  military  applications.  Air  traffic  control 
is  a  notable  civilian  application  area.  Military  applications 
range  from  air,  ocean,  and  ground  surveillance  to  missile  de¬ 
fense.  Due  to  a  drastic  increase  in  target  density  in  recent 
years,  this  subject  area  has  been  a  center  of  discussions  in 
the  open  literatures  [ 1 ]  —  [ 6 ]  and  company  reports  [7]  —  [10]  . 
Interested  readers  can  find  special  sessions  in  the  conference 
proceedings  of  recent  IEEE  Conferences  on  Decision  and  Control 
and  several  articles  for  specific  applications  in  the  IEEE 
Transactions  on  Aerospace  and  Electronic  Systems. 

Motivated  by  the  application  for  exoatmospher ic 
Ballistic  Missile  Defense  (BMD)  in  an  extremely  high  target 
density  environment,  we  have  studied  and  obtained  an  algorithm 
for  multiple  target  tracking.  References  [7]— [10]  contain 
specific  algorithms  for  exo-BMD  applications.  Due  to  the 
nature  of  this  problem  it  is  often  referred  to  as  the  scan- 
to-scan  correlation  (SSC)  problem  in  the  BMD  community. 
Although  our  algorithm  bears  this  specific  application  in 
mind,  its  concept  is  rather  general  and  it  can  be  easily  ex¬ 
tended  for  a  general  multiple  target  tracking  and  data 
correlation  application.  Our  algorithms  shares  some  features 
with  those  in  the  references,  it  also  has  some  unique 
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characteristics. 


We  will  point  out  differences  between  our  algorithm 
and  those  of  [7] -[10]  whenever  it  is  appropriate.  In  this 
section,  we  will  briefly  review  the  open  literatures  in  this 
area.  Reference  [1]  is  a  recent  and  complete  description  of  a 
multiple  target  tracking  algorithm.  Our  mild  reservations 
about  this  paper  are  (1)  it  does  not  consider  the  effect  of 
limited  sensor  resolution  and  (2)  it  attempts  to  model  every 
stage  of  the  target-tracking  process  such  as  the  a  priori 
target  distribution  and  the  probability  of  a  given  number  of 
detections  occuring  where  in  reality,  these  probabilities  may 
only  be  vaguely  known.  Reference  [2]  is  a  survey  article 
References  [3]  and  [4]  discuss  the  problem  of  track  mainte¬ 
nance  in  a  dense  target  (or  cluttered)  environment.  What  is 
missing  is  a  critical  stage  of  the  process,  track  initiation. 
The  subject  of  track  initiation  is  covered  in  [1]  and  [5]. 
Reference  [6]  discusses  a  probabilistic  data  association 
scheme  which  can  be  shown  to  be  a  special  case  of  the  al¬ 
gorithm  discussed  in  [3J— [4] .  We  will  point  out  more  spe¬ 
cifics  related  to  those  approachs  in  the  next  sections. 

This  report  is  organized  as  follows.  The  next 
section  will  give  a  very  general  discussion  of  the  hypothesis 
tree  approach  to  the  multiple  tracking  problem.  Section  3 
contains  details  of  our  algorithm.  This  algorithm  attempts  to 
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realize  the  approach  discussed  in  the  section  2  whenever  it  is 
determined  feasible.  An  example  illustrating  the  algorithm 
discussed  in  this  report  is  given  in  Section  4.  A  summary  is 
given  at  the  Section  5.  Two  appendices,  discussing  the 
polynomial  fit  formulaes  and  the  batch  estimator  used  in  this 


report,  are  given  at  the  end. 
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GENERAL  DISCUSSION 


II . 

The  multiple  target  tracking  problem  can  be  divided 
into  two  phases.  The  first  phase  is  track  initiation  and  the 
second  phase  is  track  maintenance.  They  are  discussed  indi¬ 
vidually  below. 

2 . T  Track  Initiation 

Consider  the  case  of  a  scanning  sensor.  The  first 
and  second  scan  produce  N ;  and  N2  detections,  respectively. 

The  problem  is  to  associate  the  two  sets  of  detections  to  form 
min(Ni,N/)  number  of  track  files.  Notice  that  we  have  assumed 
that  Nj  *  N;.  This  can  be  caused  by  (1)  imperfect  detection 
and  resolution,  (2)  emergence  of  new  targets  in  the  second 
scan,  and  (3)  targets  leaving  the  sensor  field  of  view  before 
the  second  scan.  In  the  following,  an  approach  for  track 
initiation  with  k  scans  of  data  is  described. 

Let  Zk  denote  all  the  measurements  (N)  collected 
during  k-th  scan,  i.e., 

zk  =  lzi(>0>  Z2(k)  , .  .  .  ,_zN(k  )  }  (2.1) 

Let  Zk  denote  the  set  of  measurements  up  to  and  including 
the  k-th  scan,  i.e., 

Zk  -  {Zi;  i=1 , • . • ,  k}  (2.2) 
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For  simplicity,  we  assume  that  N  is  the  number  of  detections 
for  all  Z^'s.  Assume  also  that  the  sensor  has  perfect  tar¬ 
get  detection.  When  this  is  not  true,  one  has  to  enumerate 
more  hypotheses  to  account  for  all  possibilities.  With  Z^, 
there  can  be  combinations  of  measurement  sequences  and 
each  measurement  sequence  represents  a  possible  track.  Let 
each  possible  combination  be  denoted  by  a  hypothesis, 

Hmk(k)  which  is  defined  by 


H 

m 


k 


(k) 


<1),z  (2) 

2 


(k)  } 


(2.3) 


Suppose  that  a  tracking  filter  is  applied  to  process  each 
possible  measurement  sequence.  The  a  posteriori  hypothesis 
probabi  lity  of  Hmk(k)  being  true  can  be  computed  recur¬ 
sively  using 


P(H  (k)/Zk) 
mk 


p(z  (k)/IJ 

nk.  rak- 1 


(k-1),Zk  ]) 


p(z  (k)/Zk"1) 
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P  (  H 
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k-1 


(k-1 ) /Zk"1 ) 


(2.4) 


k- 1 

where  p(z  (k)/H  (k-1),  Z  )  is  the  probability  density 

nk  k-1 

of  the  residual  from  the  tracking  filter  using  H  (k-1) 

mk- 1 

and  znk(k)  .  The  above  equation  can  be  derived  as  a 


*A  more  parametric  approach  for  modelling  this  probability 
density  function  is  given  in  Refs.  [1],  f 2 )  and  [5]  in  which 
situations  including  a  priori  target  distribution  and  the 
probability  of  a  given  number  of  detections  were  also  con- 
s idered , 


special  case  of  the  results  presented  in  [13],  [14].  The  final 
set  of  tracks  (total  N)  can  be  chosen  as  those  N  feasible 
hypotheses  with  the  largest  hypothesis  probabilities,  i.e., 

max  {p(Hmk(k)/Zk) ;  mk  =  (2.5) 

{N;  nmk(k  )Cf] 

where  the  feasible  set,  f  ,  is  the  restriction  that  each 
measurement  at  a  given  time  can  be  used  only  once,  i.e., 

f  ~  [Hk(k)  '  1  HjR(k)  =  $  for  i  *  j} 

The  computational  requirement  of  the  above  method  is 
clearly  non-trivial.  In  fact,  the  above  optimization 
problem  defines  a  N-d imens ional  assignment  problem.  A  well 
known  solution  to  the  2-dimensional  assignment  problem  is  the 
Hungarian  algorithm  [lb].  To  the  best  of  the  authors’ 
knowledge,  the  N-dimens ional  extension  of  the  Hungarian  al¬ 
gorithm  is  not  yet  available. 

In  many  applications,  one  may  be  able  to  pre-cluster 
the  detections  so  that  search  over  the  entire  set  of  de¬ 
tections  is  not  necessary.  Other  physical  constraints  can 
sometimes  be  imposed  to  reduce  the  search  requirements  de¬ 
pending  upon  given  systems  and  applications. 

A  similar  approach  using  a  maximum  likelihood  method 
was  described  in  [5]  in  which  the  multidimensional  search 
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problem  was  reduced  to  a  0-1  integer  programming  problem. 

2.2  Track  Continuation 

Once  track  files  have  been  established,  the  computa¬ 
tional  requirement  is  greatly  reduced.  This  is  because  for 
each  track  file  one  is  only  required  to  search  the  "admissi¬ 
ble"  region  dictated  by  the  covariance  of  the  filter  residual 
process . 

We  note  that  a  slightly  modified  method  of  the  track 
inititation  algorithm  discussed  in  2.1  can  be  applied  to  the 
track  maintenance  problem.  That  is,  one  establishes  a  new 
hypothesis  for  each  detection  resident  in  the  admissible  reg¬ 
ion.  This  procedure  results  in  an  exponentially  growing 
number  of  track  files.  One  can  inhibit  the  growing  memory  and 
computational  requirement  by  selecting  a  tree  depth  and  con¬ 
ducting  a  global  search  for  a  set  of  feasible  tracks  having 
the  highest  hypothesis  probabilities  (eqs.  (2.5),  (2.6)). 
Another  approach  is  to  combine  a  set  of  "most  likely  hypothe¬ 
ses"  growing  out  of  the  same  track  file  using  the  weighted  sum 
of  state  estimates  with  the  hypothesis  probabilities  as 
weight  ina  factors.  This  second  approach  is  the  basis  of  the 
Bayesian  tracker  presented  by  Singer  et .  al .  [3],  [4).  If  the 

depth  is  equal  to  one,  i.e.,  one  combines  all  admissible  de¬ 
tections  at  each  scan,  then  one  obtains  the  probabilistic  data 
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association  filter  of  Bar-Shalom  and  Tse  [6).  We  exmphasize 
however,  that  the  approaches  of  [3],  [4],  and  [6]  are  suitable 
for  tracking  in  a  cluttered  environment  and  do  not  directly 
address  the  multiple  target  tracking  issue. 

The  concept  discussed  above  constitutes  the  basis  of 
the  hypothesis  tree  approach  to  multiple  target  tracking.  To 
exactly  implement  the  above  algorithm  however,  will  result  in 
excessively  high  computational  requirements.  For  example, 
finding  the  optimum  solution  of  a  N-dimensional  assignment 
problem  (for  N  being  large)  is  unpractical.  A  suboptimal  but 
computationally  more  feasible  solution  is  therefore  desira¬ 
ble.  In  the  next  section,  we  present  analgorithm  which  is 
developed  specifically  for  ballistic  missile  defense  applica¬ 
tion  with  a  passive  infrared  sensor.  Several  concepts  dis¬ 
cussed  however,  are  useful  for  a  larger  class  of  multiple  tar¬ 
get  tracking  and  data  correlation  problems.  These  concepts 
will  be  identified  as  we  move  along  in  discussion. 
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ALGORITHM  DESCRIPTION 


3. 1  Introduction  and  Basic  Assumptions 


Consider  the  situation  that  there  are  N-|  and  N2 
detections  in  the  first  and  second  scans  (or  called  "frames" 
for  an  optical  sensor),  respectiely.  In  a  simple  problem  for 
which  the  target  motion  is  insignificant  between  two  scans  or 
the  relative  motion  among  targets  is  small  (such  that  the  tar¬ 
get  pattern  is  preserved),  then  one  can  apply  a  two-dimension¬ 
al  assignment  method  for  correlating  measurements  of  these  two 
scans.*  Entries  of  the  assignment  matrix  can  be  that  of 
eq.  (2.4)  with  k=2.  For  Gaussian  measurement  vectors,  one  may 
use  the  weighted  distances 


A...  =  (s.  <1)-z.  (2))T(R.(1)-rR.  (2))'1(z.  (i)-2.  (2))  (3.1) 

1  J  J  *•  J  1  J 

as  entries  where  _Z£(k)  is  the  i-th  measurement  of  the  k-th 
scan  with  measurement  covariance  R^(k).  Once  measurement  of 
two  scans  have  been  correlated,  a  velocity  vector  can  be 
established  making  the  correlation  with  measurements  of  the 
third  frame  somewhat  easier.  For  an  optical  sensor  measuring 
1  ine-of-s ight  angles,  this  velocity  vector  will  initially  be 
limited  to  the  angle  domain  while  a  radar  sensor  can  give  a 
three  dimensional  velocity  vector. 

The  tracking  pioblem  discussed  in  this  report  is 


*This  is  similar  to  a  two  sensor  measurement  correlation 
problem  discussed  in  [11]. 
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more  complicated  than  that  described  above.  Target  motion, 
density  and  limited  sensor  resolution  are  such  that  target 
patterns  are  not  preserved  in  successive  scans.  In  this  case, 
one  must  apply  more  knowledge  about  the  target  motion  dynamics 
and  use  more  scans  of  data  to  identify  a  string  of  successive 
measurements  representing  the  same  target. 

In  the  following  two  subsections,  we  will  discuss 
the  problem  of  track  initiation  and  continuation  individu¬ 
ally.  An  overall  description  is  given  in  Table  3.1. 

3.2  Track  Initiation 

The  most  crucial  and  difficult  part  of  the  multiple 
target  tracking  problem  is  track  initiation.  Specifically  for 
the  optical  sensor  tracking  problem,  the  following  factors 
further  complicate  the  issues: 

(1)  The  target  angular  velocities  vary  over  a  wide 
range  of  values  making  the  acceptance  cell  on  the 
second  frame  large  resulting  in  a  large  number  of 
false  correlations. 

(2)  With  angle-only  (passive  receiver)  measurements, 
the  target  range  estimate  can  not  be  readily  ob¬ 
tained  making  impossible  the  use  of  precise 
ballistic  equations  of  motion  as  target  dynamics. 

In  the  case  of  radar  tracking,  the  first  point  above 

may  still  be  true  depending  on  target  velocity  and  data  rate. 

The  second  point  above  is  at  least  partially  true  since  at  the 

initial  stage,  a  large  number  of  track  files  are  false,  using 
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TABLE  3.1 


ALGORITHM  DESCRIPTION 


junctions  Description 

Initiation- 


No.  of  Scans 

First  Frame  Angular 
Velocity 

Prediction 

Correlation 

Continuation 


5-8 

A  Priori  Knowledge  for  some 
targets,  recursively  search 
for  parallel  targets. 

Oth  -  3rd  order  polynomial 

(1)  Track  split 

(2)  Chi-square  test 

(3)  Pattern  match  test 


Prediction 
Search  Bin  Size 
Correlation 


Target  equation  of  motion  or 
polynomial  dynamics 

Filter  covariance  and  model 
error  analysis 

(1)  Track  split 

DoHAl*rt  1  .  .  J 

\w  i  uuci  ii  main i  itJM 
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the  exact  target  dynamics  at  this  time  is  very  time  consuming. 
For  these  reasons,  we  use 

(1)  a  general/parallel  search  scheme  for  reducing  the 
number  of  false  correlations  and 

(2)  a  polynomial  function  in  each  angle  domain  as  target 
dynamics  to  simplify  calculations. 

Furthermore,  a  sliding  window  scheme  is  employed  to  initiate 

tracks  for  new  detections  in  each  frame.  It  is  important 

to  note  that  this  is  an  iterative  rather  than  recursive  method 


and  data  from  many  frames  must  be  saved. 

The  track  initiation  logic  is  illustrated  in  Fig. 
3.1.  We  use  the  following  steps  to  illustrate  the  general  and 
parallel  search  scheme. 


(1  ) 


For  a  given  detection  in  the  first  frame  (called  an 
initiator),  draw  an  acceptance  region  centered  at 
this  detection.  Detections  on  the  second  frame 


which  fall  into 
The  size  of  the 


this  region  form  potential  tracks, 
initial  acceptance  region  is  deter— 


mined  by  the  maximum  target  angular  velocity. 
Usually  a  large  number  of  potential  tracks  result 
for  one  given  initiator. 


(2)  Apply  the  straight  line  extrapolation  scheme  (see 

Appendix  A)  to  extend  all  potential  tracks  into  the 
third  frame.  The  size  of  the  acceptance  region  at 
the  third  frame  is  determined  by  model  and  measure¬ 
ment  errors  of  the  linear  extrapolation.  This 
acceptance  region  size  is  usually  much  smaller  than 
that  of  the  first  step  above. 


(3)  Apply  a  second  order  polynomial  (see  Appendix  A)  to 
extend  tracks  into  the  fourth  frame.  Similarly,  the 
acceptance  region  size  further  reduces.  Tracks 
which  do  not  receive  a  detection  in  their  acceptance 
region  will  be  dropped.  Tracks  receiving  multiple 
measurements  in  their  acceptance  region  will  be 
split. 
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Fig.  3.1.  Track  initiation  logic. 
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(4)  Continue  to  a  total  of  5  to  8  frames  depending  upon 
the  target  density  and  scenario.  At  the  end  of  this 
stage,  usually  only  a  few  potential  tracks  remain. 
Final  choices  of  tracks  are  selected  using  a  global 
polynomial  fit.  Those  tracks  witn  weighted  residu¬ 
als  (Chi-squares)  below  a  threshold  will  all  be  re¬ 
tained.  This  completes  the  general  search  for  an 
initiator. 

(5)  Go  back  to  the  first  frame,  assuming  that  targets  in 
the  neighborhood  of  the  initiator  will  travel  in 
nearly  the  same  direction;  one  therefore  only  has  to 
search  for  detections  in  the  successive  frames  in 
parallel  with  the  track(s)  established  with  this 
initiator.  This  step  greatly  reduces  the  computa¬ 
tional  and  memory  requirements.  This  step  is  called 
parallel  search. 

(6)  Once  all  parallel  tracks  have  been  found,  go  back  to 
the  first  frame,  find  another  detection  which  has 
not  been  included  in  any  tracks  to  use  as  a  new 
initiator  for  general  search. 

(7)  Repeat  until  all  detections  of  the  first  frame  have 
been  exhausted. 

Once  the  initial  correlation  described  above  has 
been  completed,  one  can  trim  track  files  by  applying  a  n-th 
order  polynomial  (n  is  determined  by  a  particular  application) 
fit  to  measurements  of  a  file  and  reject  those  files  with 
excess 'vely  high  ’'chi-square"  values.  A  third  correlation 
method  listed  in  Table  3.1  under  track  initiation  is  called 
"pattern  match  test".  This  test  is  the  same  as  the  one  used 
in  the  track  continuation.  We  therefore  defer  its  explanation 
until  the  next  subsection. 

The  above  procedure  is  applied  in  a  sliding  window 
fashion  so  that  measurements  not  included  in  the  track  file 
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are  used  to  initiate  new  tracks.  This  is  illustated  in  Fig. 
3.2. 

3.3  Track  Continuation 

The  track  initiation  process  correlates  measurements 
over  5  to  8  frames  to  produce  track  files  using  data  of  all 
frames  (therefore  a  smoothing  process).  The  track  continua¬ 
tion  stage  can  be  a  traditional  prediction,  correlation  and 
updating  process.  As  shown  in  Table  3.1,  the  prediction  step 
may  use  either  target  equations  of  motion  or  polynomial 
dynamics  for  reducing  the  computational  burden.  In  the 
exoatmospheric  BMD  application,  we  have  found  that  the  use  of 
a  precision  tracking  filter  with  the  complete  target  equations 
of  motion  greatly  enhances  the  performance.  We  have 
implemented  both  the  extended  Kalman  filter  and  the  batch 
filter  described  in  [16]  for  track  continuation .  The  batch 
filter  is  also  briefly  reviewed  in  Appendix  B  for  the  purpose 
of  quick  reference. 

We  use  Figure  3.3  to  illustrate  some  typical  situa¬ 
tions  encountered  in  the  track  continuation  process.  As  a 
track  moves  along,  if  multiple  detections  are  encounted,  the 
track  is  split  (case  i).  If  no  detections  are  found  for 
several  frames  in  a  row,  the  track  is  dropped  (cases  1  and 
2).  There  may  also  be  situations  for  which  a  track  is  split 
and  then  merged  (case  3). 
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Fig.  3.2.  Sliding  window  initiation  logic. 
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Ambiguities  may  arise  when  a  number  of  track  files 
have  overlapping  acceptance  regions  and  share  the  same 
detections.  This  is  illustrated  in  Figure  3.4  for  the  case 
where  two  tracks  share  two  measurements.  Problems  of  this 
kind  are  similar  to  the  assignment  problem  in  operations 
research.  The  optimal  (in  the  sense  of  minimum  sum  of 
weighted  distances)  solution  is  usually  obtained  with  a 
so-called  Munkres'  algorithm  (see  [11]  and  [15]).  If  one 
attempts  to  resolve  this  ambiguity  at  the  frame  where  it  is 
encountered,  this  is  called  immediate  conflict  resolution. 
Since  the  track  formation  problem  is  really  a  multiple 
dimensional  assignment  problem  (eq.  (2.5)),  a  more  reliable 
decision  can  be  obtained  by  deferring  decisions  until  further 
measurements  have  been  received.  This  is  called  deferred 
conf 1  ict  resolution.  A  tradeoff  for  these  methods  is 
computation  vs  performance.  We  have  implemented  both 
the  immediate  conflict  resolution  method  and  the  one  frame 
deferred  conflict  resolution  method.  Later  numerical  results 
will  compare  the  performance  of  these  two  methods. 

We  use  Fig.  3.5  to  further  illustrate  the  conflict 
resolution  methods.  In  Fig.  3.5a  two  track  files  with 
measurements  up  to  the  N-1  frame  are  extrapolated  to  Nth  frame 
and  found  to  compete  for  measurements  a  and  b.  To  resolve 
this  conflict  immediately  is  to  first  form  a  distance  matrix 


A  Immediate  Resolution 


Frame#  N 
A  Delayed  Resolution 


N  +  2 


+  Using  Interpolated  Estimate  to  Resolve  the 
Ambiguity  of  the  Previous  Frame 


Fig.  3.4.  Track  continuation  ambiguity  resolution. 


TR-643( 3. 5a) 


Track  #  1 


2 


a  and  b  in  common 
admissible  region 


/  V 

s  \ 

/  \ 


\ 


\ 


b 


Frame  #  N-l  N 


Pattern  Match  Matrix  • 

Measurements 
a  b 

=*=  1 

O 
<T3 

L. 

*-  2 


Fig.  3.5(a).  Illustration  of  ambiguity  resolution  technique: 
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la 


D 


2a 


D 


D 


lb 


2b 


20 


|T*-643(3.5b) 


Track  #  1 


2 


Track  Split 
at  N-th  Frame 


c  and  d  in  common 

admissible  region 

— — • 

✓ 

/  c 

/ 


Frame  #  N-l  N  N  +  l 


Two-Layered  Pattern  Match  Matrix  : 


For  Measurement  c 
Measurements 


For  Measurement  d 
Measurements 


a 

b 

a 

b 

1 

Dlac 

Dlbc 

D  lad 

Dlbd 

2 

D2ac 

D2bc 

« 

°2ad 

°2bd 

■ 

Fig.  3.5(b).  Illustration  of  ambiguity  resolution  technique: 
multiple-layered  pattern  matching. 
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Fig.  3.5(c).  Illustration  of  ambiguity  resolution  technique; 
an  one-frame  deferred  resolution  method. 


22 


shown  in  Fig.  3.5a.  For  example,  D-]a  is  the  weighted 
distance  between  the  predicted  location  of  the  first  track  and 
measurement  a.  The  "optimum"  solution  is  the  pair  of  D^j's 
that  gives  the  minimum  sum  {see  for  example,  [11]).  This  is 
the  "immediate"  resolution  method. 

If  one  simply  splits  all  tracks  and  moves  to  the 
N+lst  frame  (Fig.  3.5b)  and  assumes  that  there  are  again  two 
measurements  (c  and  d)  falling  into  the  common  admissible 
region,  there  can  be  a  total  of  eight  possible  tracks  formed. 
The  "optimum"  approach  for  this  problem  is  to  first  form  a 
"two-layered"  matrix  with  dimension  (2x2x2).  The  entries  of 
this  matrix  are  the  posteriori  hypothesis  probabilities  shown 
in  (2.4).  The  "optimum"  solution  is  the  pair  of  entires  with 
minimum  sum.  Two  factors  will  have  to  be  considered:  (1)  the 
computation  of  (2.4)  is  tedious  especially  when  the  numbers  of 
track  files  and  measurements  involved  are  large  and  (2)  an 
efficient  algorithm  (equivalent  to  the  Hungarian  Method  for 
searching  the  optimum  solution)  has  not  been  found.  For  the 
above  reasons,  we  device  a  suboptimal  procedure  which  is 
particularly  efficient  when  the  numbers  of  track  files  and 
measurements  are  large.  Notice  that  there  are  two  possible 
tracks  to  go  from  track  file  #1  through  measurement  a  of  the 
Nth  frame,  namely  (1,a,c)  and  (1,a,d)  (see  Fig.  3.5c).  We 
first  select  a  track  among  these  two  with  the  smallest  residu- 
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al .  This  procedure  is  repeated  for  all  track  files  and  all 
measurements  in  the  Nth  frame.  The  residuals  at  the  Nth  frame 
(a  smoothed  residual)  of  all  selected  tracks  are  used  to  form 
the  distance  matrix  (Fig.  3.5c).  The  final  selection  is  the 
pair  of  tracks  with  the  minimum  sum  of  entries  of  the  distance 
matrix.  This  is  the  one-frame  deferred  resolution  method. 
Notice  that  the  one  frame  deferred  method  is  analogous  to  the 
one-step  lagged  fixed-lag  smoothing. 

Notice  that  the  above  is  attempting  to  reduce  a 
N--d imensional  (N--3  for  this  case)  assignment  problem  to  a 
2-dimensional  assignment  problem.  This  procedure  although 
suboptimal,  is  straightforward  to  implement  and  gives  better 
performance  than  the  immediate  resolution  method  (see  section 
4)  . 

The  deferred  ambiguity  resolution  method  can  be  imple¬ 
mented  in  a  sliding  window  fashion,  i.e.,  one  uses  data  in  the 
past  and  future  for  resolving  the  conflict  of  the  current 
frame.  This  method  is  also  used  by  the  track  initiation  pro¬ 
cess  right  after  the  "Chi-square"  test  (see  Section  3.2  and 
Table  3.1).  It  is  referred  to  as  the  pattern  match  test  in 
Table  3,1. 

3.4  Track  File  Smoothing 

The  track  initiation  algorithm  described  above  is  a 
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smoothing  algorithm  because  it  uses  all  past  data  in  a  batch 
processing  mode.  All  past  data  will  therefore  have  to  be 
stored  for  this  purpose.  If  one  exactly  implements  eq,  (2.4) 
for  track  initiation,  this  process  can  be  recursive  although 
the  number  of  tracks  can  grow  out  of  hand  rather  quickly.  The 
track  continuation  algorithm  using  the  one  frame  deferred 
resolution  method  is  a  one  step  lagged  fixed-lag  smoothing 
process;  one  is  only  required  to  store  the  immediate  past  set 
of  measurements. 

In  the  exoatmospher ic  BMD  application  where  a 
passive  optical  (Infrared)  sensor  is  used  for  target 
tracking, a  conventional  recursive  filter  (e.g.,  the  extended 
Kalman  filter)  may  not  work  satisfactorily,  [16].  Instead,  a 
maximum  likelihood  estimator  based  batch  estimator  is  found  to 


give  near  optimum  performance  (see  [16]  and  also  Appendix  B  of 
this  report  for  the  batch  algorithm) .  This  estimator  however 
requires  that  all  measurements  be  saved.  Suppose  that  all 
past  measurements  have  been  saved,  then  one  can  further  edit 
track  file  measurements  when  the  batch  filter  is  being  used  to 
process  these  track  files.  This  process  is  illustrated  in 


Pig.  3.6.  A  hypothesis  test  is  applied  to  measurements  of  a 


track  file  with  respect  to  the  estimated  trajectory.  When  a 


residual  is  too  large,  that  particular  measurement  is  rejected 


and  a  new  measurement  is  chosen.  This  procedure  indeed  works 
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well  as  will  be  illustrated  in  the  numerical  results  section. 


3.5  A  Partitioned  Field  of  View  Implementation 

From  the  discussion  of  previous  sections,  it  is 
clear  that  the  solution  to  the  multiple  target  tracking 
problem  is  not  a  problem  requiring  sophisticated  mathematical 
manipulations,  rather  involving  simple  calculations  and  large 
scale  data  and  file  management  schemes.  When  it  is  to  be 
implemented  on  a  digital  computer,  the  problem  of  indexing 
measurements  in  a  track  file  can  be  a  real  challenge  for 
programmers . 

In  this  section,  we  briefly  describe  an 
implementation  scheme  which  is  found  to  be  computationally 
efficient.  The  sensor  field  of  view  (FOV)  is  first 
partitioned  into  bins  (Figure  3.7).  In  the  exoatmospheric  BMD 
appx icat ion ,  these  will  specifically  be  bins  in  azimuth  and 
elevation.  The  bin  size  should  not  be  smaller  than  the 
maximum  target  motion  between  two  observation  and  should  be  at 
least  a  few  standard  deviations  of  prediction  error.  When 
measurements  of  a  bin  are  being  track  initiated,  only 
measurements  in  the  bin  itself  and  the  bordering  bins  for  the 
subsequent  measurement  frame  are  searched.  In  the  track 
continuation  mode,  only  those  track  files  with  their  last 
measurement  residing  in  the  current  bin  and  its  immediate 
neighbors  are  used  for  processing.  For  the  target  density 
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Fig.  3.7.  Partitioned  field  of  view  implementation. 
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considered  in  the  BMD  application,  this  method  significantly 
reduces  memory  access  time.  Furthermore,  if  one  stores  all 
measurements  and  track  files  on  disk,  then  the  core  memory 
will  only  have  to  be  large  enough  to  hold  those  of  a  few  urns. 

Another  advantage  of  the  partitioned  field  of  view 
approach  is  that  it  is  especially  suitable  for  implementation 
with  dedicated  multiple  parallel  processors.  In  this  case, 
each  bin  may  be  handled  by  an  independent  processor  and  a 
monitor/supervisor  processor  can  be  assigned  to  handle  tracks 
crossing  bin  boundaries. 

3.6  Summary 

In  this  section,  we  have  described  a  multiple  target 
tracking  algorithm  with  specific  applications  to  ballistic 
missile  defense.  Features  of  this  algorithm  are  summarized 
below. 

(1)  The  track  initiation  process  uses  a  general/parallel 
search  procedure  which  can  substantially  reduce  pro¬ 
cessing  time. 

(2)  The  track  initiation  applied  to  the  exoatmospheric 
BMD  problem  requires  5-7  frames  of  measurements.  It 
is  applied  in  a  sliding  window  fashion  to  handle  the 
changing  scene  and  crossing  traffic  problem, 

(3)  A  one  frame  deferred  conflict  resolution  scheme 

is  presented  for  resolving  the  problem  of  multiple 
track  files  competing  for  several  measurements. 

(4)  If  all  or  part  of  the  past  measurements  are  being 
saved,  a  concept  using  a  precision  filter  for 
further  track  file  editing  is  described. 
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(5)  A  partitioned  field  of  view  implementation  scheme, 
with  potential  for  multiple  processor  implementa¬ 
tion,  is  discussed. 
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IV. 


EXAMPLE 


In  this  section,  we  illustrate  our  results  using  an 
example.  This  example  represents  a  typical  high  target  densi¬ 
ty  environment  in  BMD  systems. 

The  exoatmospher ic  BMD  system  concept  calls  for  an 
infrared  sensor  deployed  on  a  probe  vehicle  flying  on  a  bal¬ 
listic  trajectory  to  observe  and  track  incoming  targets.  In 
the  example  presented  in  this  section,  a  complex  consisting  of 
approximately  600  targets  are  being  observed  with  20  frames  of 
measurements  with  ten  seconds  of  time  between  frames. 

The  target  true  angle  data  at  each  frame  are  first 
processed  by  a  sensor/signal  processor  (SSP)  simulation  to 
generate  simulated  angle  measurements.  This  simulation  takes 
into  account  1)  target  intensity  variations,  2)  background  and 
receiver  interference,  and  3)  effects  due  to  limited  sensor 
resolution.  Targets  may  therefore  be  detected  in  one  frame 
but  missed  in  the  next.  Several  closely  spaced  targets  may 
not  be  mutually  resolved  and  therefore  result  in  fewer  number 
of  measurements  than  targets.  Measurements  containing  unre¬ 
solved  targets  are  usually  biased  with  a  standard  deviation 
determined  by  all  targets  involved.  Some  resolved  targets  may 
also  contain  excessively  high  measurement  errors  due  to  inter¬ 
ference  introduced  by  nearby  targets.  A  detailed  description 
of  these  effects  can  be  found  in  [17]. 
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Tracking  performance  is  evaluated  using  a  track  file 
consistency  measure-  We  first  illustrate  some  typical  track 
files  using  Table  4.1.  The  top  row  gives  frame  number.  The 
left  column  gives  track  file  (TF)  identification  number.  The 
entries  are  target  identification  numbers.  Only  12  frames  of 
data  are  shown  for  illustration.  Track  file  (TF)  #100  con¬ 
tains  target  20  throughout  12  frames  shown.  This  is  a  well- 
defined,  consistant  track  file.  TF  #101  contains  target  22. 

It  begins  at  frame  4  and  misses  the  target  at  frame  7.  TF 
#102  and  103  should  be  examined  together.  Targets  31  and  32 
form  an  unresolved  closely  spaced  target  cluster  at  frame  1, 

2,  3,  4,  and  7.  When  these  two  targets  begin  to  get  resolved 
at  frame  5,  6,  and  beyond  frame  7,  these  two  track  files  have 
successfully  tracked  them.  TF  #104  and  105  show  a  similar 
situation  except  that  they  have  failed  to  consistently  track 
the  CSO  splitting  (see  frame  8  and  beyond).  Unresolved 
measurements  containing  more  than  two  targets  can  also  be 
found  in  the  threat  data  examined. 

Based  upon  the  above  illustration,  we  now  describe  a 
"target  oriented"  scoring  (performance  evaluation)  scheme. 

For  a  given  target,  we  first  identify  all  track  files  contain¬ 
ing  it.  For  example,  using  Table  4.1  for  target  number  31, 
one  finds  that  both  track  files  102  and  103  contain  this  tar¬ 
get.  A  track  file  containing  this  target  most  often  is 
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TABLE  4.1 


ILLUSTRATION  OF  TYPICAL  TRACK  FILES 


'^s\Frame  # 

TF 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

100 

20 

20 

20 

20 

20 

20 

20 

20 

20 

20 

20 

20 

101 

0 

0 

0 

22 

22 

22 

0 

22 

22 

22 

22 

22 

102 

31 

31 

31 

31 

31 

31 

31 

31 

31 

31 

31 

31 

32 

32 

32 

32 

32 

103 

31 

31 

31 

31 

32 

32 

31 

32 

32 

32 

32 

32 

32 

32 

32 

32 

32 

104 

26 

26 

26 

26 

26 

26 

26 

27 

27 

27 

27 

27 

27 

27 

27 

27 

27 

105 

26 

26 

26 

26 

27 

27 

26 

26 

26  . 

26 

26 

26 

27 

27 

27 

27 

27 

Entries  are  target  ID  s. 
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assigned  to  represent  this  target.  In  Table  4.1,  track  files 
100,  101,  102,  103,  104,  and  105  are  assigned  to  represent 
targets  20,  22,  31,  32,  27,  and  26,  respectively.  It  is  clear 
that  a  track  file  can  sometimes  be  assigned  to  more  than  one 
targets  because  of  the  closely  spaced  target  effect.  A  per¬ 
formance  for  each  "target"  can  then  be  evaluated  using  the 
assigned  track  files.  Suppose  that  the  number  of  true  occur¬ 
rences  of  each  target  in  Table  4.1  is  the  same  as  the  table 
shows,  then  we  conclude  that  targets  20,  22,  31,  32,  26,  and 
27  each  meet  100%,  100%,  100%,  100%,  83.33%  and  83.33%  per¬ 
formance,  respectively.  Take  target  26  as  an  example.  Target 
26  has  appeared  for  12  frames.  The  track  file  assigned  to 
represent  target  26  is  track  file  number  105  which  contains 
target  26  for  10  frames.  This  means  that  the  performance  on 
tracking  target  26  is  10/12  =  83.33%. 

We  apply  the  above  scoring  scheme  to  tracks  genera¬ 
ted  for  the  600  target  case  described  earlier,  the  results  are 
shown  in  Fig.  4.1.  The  horizontal  axis  gives  performance  cri¬ 
terion  described  in  the  previous  paragraph  and  the  vertical 
axis  gives  the  percentage  of  targets  satisfi.ng  a  given  crite¬ 
rion.  Again  using  data  of  Table  4,1  for  the  purpose  of  illus¬ 
tration,  4  out  of  6  targets  meet  the  100%  performance  criteri¬ 
on  and  all  targets  meet  the  80%  criterion.  We  use  results 
shown  in  Fig.  4.1  to  compare  the  three  tracking  algorithms  de- 
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Fig.  4.1.  Illustration  of  tracking  performance. 
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scribed  in  the  previous  section.  Clearly,  the  most 
sophisticated  algorithm  gives  the  best  result.  Also  notice 
that  the  biggest  gain  in  using  a  sophisticated  algorithm  is  at 
the  100%  performance  criterion  level  while  the  differences  at 
the  lower  criterion  region  grow  smaller.  This  is  because  the 
precision  trajectory  estimation  begins  by  using  track  files 
produced  by  the  correlation  algorithm.  In  order  for  the 
precision  estimation  to  succeed,  track  files  presented  by  the 
correlation  algorithm  will  have  to  be  reasonably 
cons  is  tent ( e . g .  ,  >  50%).  The  gain  shown  in  Fig.  4.1  for  the 
case  with  track  file  smoothng  is  obtained  largely  by  upgrading 
the  track  files  satisfying  60%  -  70%  performance  criteria. 

In  Table  4.2,  we  compare  the  processing  time  of 
these  three  methods.  Notice  that  they  are  compared  on  a  rela¬ 
tive  basis  and  the  total  time  includes  both  correlation  and 
precision  tracking. 
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TABLE  4.2 

PROCESSING  TIME  COMPARISON 


CPU  Time 

(normalized) 

Immediate  Resolution 

1 

Delayed  Resolution 

1.2 

With  Track  Fiie 

1.8 

Smoothing 


Notes:  Processing  time  includes  correlation  and 
precision  estimation. 
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V. 


SUMMARY 


In  this  report,  we  have  described  an  algorithm  for 
multiple  target  tracking  and  data  correlation  in  a  dense  tar¬ 
get  environment.  Although  some  of  our  discussions  were  cen¬ 
tered  around  the  tracking  problem  with  a  passive  infrared  sen¬ 
sor,  the  approach  represented  by  this  report  is  applicable  in 
a  general  situation. 

The  significant  features  of  our  approach  include: 

1)  the  use  of  a  general/parallel  search  for  track  initiation 
and  2)  the  one  frame  deferred  ambiguity  resolution  method.  A 
partitioned  field  of  view  implementation  scheme,  with  poten¬ 
tial  for  multiple  processor  implementation,  was  also  discus¬ 
sed  . 


38 


APPENDIX  A:  Some  Useful  Results  in  Fitting  Polynomials  to  a 
Set  of  Noisy  Measurements 

The  techniques  of  polynomial-fit  have  found  wide 
acceptance  in  applications.  In  this  Appendix,  we  present  a 
brief  analysis  of  some  polynomial-fit  formulas  which  are 
relevent  to  trajectory  tracking  and  prediction  problems.  The 
polynomial  predictor  used  in  the  scan-to-scan  correlation 
algorithm  can  be  derived  using  these  results. 

A. 1  General  Results 


Consider  a  set  of  polynomials  of  variable  t 
(denoting  time  in  our  applications)  represented  by  f0(t), 
f  1  ( t fjc(  t )  .  Assume  that  this  set  of  polynomials  can 
adequately  represent  a  function  y(t)  in  the  form  of  a  weighted 
linear  sum,  i.e.. 


y(t) 


akfk(t). 


(A.  1  ) 


Given  a  set  of  noise  corrupted  discrete  time  measurements  of 
y(t) 

z(tn)-y(tn)  +  V  <A-2) 

where  5n  is  an  uncorrelated  non-stat ionarv  zero  mean  noise 
sequence  with  variance  on2.  The  objective  is  to  find  a 
set  of  a^ ' s  giving  the  "optimal"  estimate  of  y(t)  from 
z(tn),  n=1,...,N.  If  the  weighted-least-square  estimator  is 
used,  then  the  optimal  aj^'s  are  those  minimizing 


39 


i 


N  1 


(z(tn)  -  y(tn) ) 


(A. 3) 


Let 


a  •- 


L“k 


and  £  = 


Lf„J 


(A. 4) 


Then  the  estimate  of  a,  is 
N 


n-1  n 
The  covariance  of  a,  Pf  is 

N 


-1  .  N 

i-  fE  -T3  i(tn>  £t“„>]  [E  "V  ' 

L- n  =  1  n 


P  = 


E  71  i(t")  lT(t": 

n=1  n 


-1 


(A. 6) 


Consider  a  special  set  of  polynomials 

fQ(t)  =  1 
f  1  <  t )  =  t 


(A. 7) 


fR(t)  =  tK/K! 


Then , 
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F(tn)  - 


tnK/K! 


*  tn  .  tnVKl 


tn . tnK/K! 


tnK/KI  .  tn2K/(KI)2 


(A. 8) 


and 


i-lfiW  Bi: 

-IS  ~ 


2  Z(tn)  £<fcn> 


2  F(V 


(A. 9) 


(A. 10) 


A. 2  Results  for  a  Second  Order  Polynomial  With  Uniform  Time 
Samples  and  Stationary  Noise  Sequences 


Consider  the  following  conditions 
'•  fK(t)  *  £2(t)  *  1 

2.  Time  samples  are  taken  uniformly  and  t=Q 

corresponds  to  the  center  of  the  data  batch 
This  implies 

( N- 1 ) T 


t  =  ( n- 1 )T  -  „ 

n  2 


(A, 11) 
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where  T  is  the  intersample  spacing. 


3.  0n  =  o  for  all  n's 


Also,  using  the  following  identities, 


N 


1.  £  <  ( n  —  1  )  T  - 

n  =  1 


( N- 1 )  T 


-  )  =  0 


N 


L  (  (n-DT  -  — — )  =  j~2  (N-1)N(N+1) 
n  =  1  "  ' 


2  2 


(A. 12) 


N 


3.  £  ( (n-1 )T  - 

n  =  1 


( N-  1  )  T 
2 


3 

)  =  0 


E.  ((n-.)T  -  )4  -  Sli£=m<NlU  ,3N2-7) 


n=  1 


Then  Equations  (A. 9)  and  (A. 8)  become 


P=o“ 


3  _ (  3N  -7 ) 

4  (N-2)N(N+2) 


0 

12 


-30 


-30 


T  (N-1 )N(N+1 ) 


T  (N-1 )N(N+ 1 ) 
0 


T  ( N-2 ) N( N+2 ) 
0 


720 


T  (N-2) (N-1 )N(N+1 ) (N+2) 


and 


(A. 13) 
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N 


N 

s, 


”1 


z(t  ) 
n ' 


T(  n- 1  ^2-~)  z(tn) 


T2(n_1„iJi_lI)  z(tn)/2 


(A. 14) 


If  one  repeats  the  above  exercise  by  choosing  K=1,  one  will 
obtain 


P  =  o 


2 

N 


0 


12 


T2{N-1 }N(N+1 ) 


(A. 15) 


i, 


z(t 


N 


V  T(  n  - 1  _ 

2  ]  Z(tr 


(A. 16) 


A. 3  Applications  to  Position,  Velocity  and  Accelera¬ 
tion  Estimation 


Let  p0,v0,  and  a0  denote  the  position,  velocity, 
and  acceleration,  respectively,  of  a  moving  object  at  t=0, 


43 


then  its  position  at  an  arbitrary  time  t  is 

P(t)  =  pQ  +  vQt  +  aQt2  (A.  17) 

If  p(tn),  n=1,...,  N  denote  a  set  of  noisy  measurements  of 
p(tn),  the  Equations  (A. 13),  (A. 14),  (A. 15),  and  (A. 16)  can 
be  directly  applied  for  estimating  p0,v0,  and  aQ. 

Equation  (A. 13)  and  (A. 14)  correspond  to  a  constant  accelera¬ 
tion  model  while  Equations  (A. 15)  and  (A- 16)  corresponds  to  a 
constant  velocity  model.  With  p0,v0,  and  a0,  estimates 
at  time  t  are  obtained  by  using 

-  P(t)-|  r  PQ 

v(t)  =  4  (t)  vQ 

_  a{ t )  L  ao 

where 

1  t  t2/2 

0  1  t  (A. 19) 

0  0  1 

and  the  covariance  becomes 

P(t)  =  4>(t)  P  $T(t)  (A. 20) 

The  results  for  a  constant  velocity  model  can  be 
obtained  accordingly. 
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APPENDIX  B:  Cn  the  Multiple  and  Single  Stage  Iterative 
Least  Square  Estimators 

The  batch  state  estimator  based  cn  the  Maximum  Like¬ 
lihood  or  the  Weighted  Minimum  Mean  Square  error  criterion  for 
the  angle  only  tracking  application  was  discussed  in  detail  in 
[16].  In  this  appendix,  we  briefly  state  this  algorithm  in 
section  B. 1 .  It  is  called  a  multiple  stage  iterative  al¬ 
gorithm  because  it  attempts  to  minimize  residuals  with  respect 
to  several  (all)  measurements.  If  one  selects  to  minimize 
with  respect  to  the  most  recent  measurement  only  while  holding 
the  estimate  obtained  with  all  previous  measurements  constant, 
it  is  called  the  single  stage  iterative  algorithm.  This  al¬ 
gorithm  is  discussed  in  section  B.2.  Several  forms  of  itera¬ 
tive  filters  can  be  found  in  [18]. 

B . 1  A  Multiple  Stage  Iterative  Least  Square  Estimator 
Consider  the  following  state  and  measurement  equa- 

t ions  : 

x  -  f(x)  ( B. 1 ) 

Zy^  =  h(x^)  +  y_k;  k=1,...,N  (B.2) 

where  x  is  the  state  vector,  zj^  is  the  measurement  vector, 
and  vk  is  the  measurement  noise  vector  with  Gaussian  dis¬ 
tribution  with  zero  mean  and  known  covariance.  The  current 
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time  is  denoted  by  the  index  k=N.  We  state  the  estimator 
equations  without  derivation. 


Let  x„  denote  the  k-th  iteration  of  the  estimate 
—N/N 


of  xN,  then 


' k+ 1  ~k  k 

-N/N  ~  -N/N  *N/N 


where 


n 

rk 

GN 

,k 

*n 

1  s 


-1 


’N/N 


E 

n=  1 
N 

-  E 


n=1 


T  T 

-k  uk  _-1.  "k  ,, 

Gn  Hn  Rn  — n  '  ^Vn11 


Gk  Hk  R 
N  n  n  n  n 


k'1  k 

Gn+1  '  n»N-1,  N-2 


I  (an  identity  matrix) 
linearized  transition  matrix 


Lne  Jacobian  matrix  of 


(B.3) 


(B.4) 


n 

'k 

-n/N 


the  Jacobian  matrix  of  h(*n/N) 


result  of  integrating  backward  from 

tN  tQ  fcn  USin9  *N  =  —N/N " 


'k+1  ~k 

The  iteration  is  terminated  when  xvl and  x.,  are  suffi- 

— N/N  —N/N 

ently  close.  The  notation  Xjyj  denotes  the  estimate  of 


based  upon  data  z 


i 


We  make  the  following  remarks: 


1)  The  above  algorithm  is  a  realization  of  the  maxi¬ 
mum  likelihood  estimator  with  Gaussian  measurement 
noise  process.  It  is  well-known  that  the  maximum 
likelihood  estimate  is  asymptotically  efficient 
and  Gaussian  and  approaches  the  Cramer-Rao  bound. 


2) 


The  Pn/jsj  of  (B.4)  is  an  approximate  expression 


for  the  covariance  of  x_N^N. 


The  PN/N  evaluated 


at  the  true  state  is  the  Cramer-Rao  lower  bound  on 
the  covariance  of  Since  2Ln/n  approaches 

the  true  state  with  probability  one,  also 

approaches  the  Cramer-Rao  bound  with  probability 
one . 


3)  Notice  that  the  inverse  of  P>j/n  is  the  Fisher's 
information  matrix.  The  invertibil ity  of  the 
information  matrix  is  tied  to  the  observabilil ity 
of  the  system,  see  for  example  [19], 


There  are  many  application  areas  for  this  algorithm.  One 
important  application  area  is  track  initiation.  Since  the 
initial  covariance  and  state  estimates  are  not  generally 
given  a  priori,  the  above  algorithm  can  obtain  the  best 
estimates  based  on  the  first  N  measurement  vectors  and 
then  proceed  to  use  x^/n  anc^  PN/N  as  initial  state  and 

covariance  estimates,  respectively.  This  method  is  some¬ 
times  referred  to  as  the  information  matrix  approach  for 
filter  initiation. 


47 


B.2  A  Single  Stage  Iterative  Least  Square  Estimator 


The  above  results  were  given  without  derivations. 
Interested  readers  can  consult  [16]  for  details.  Notice  that 
the  estimator  (eqs.  (B.3))  is  a  weighted  combination  of  all 
previous  measurements.  This  requires  storing  of  all  past 
measurements.  The  well-known  extended  Kalman  filter  only  re¬ 
quires  storing  the  most  recent  measurements.  However  it  may 
result  in  biased  estimates.  If  indeed  past  measurements  can 
not  be  stored,  one  can  extend  the  result  of  Section  B. 1  to  a 
single  stage  iterative  estimator.  It  can  also  be  shown  that 
the  extended  Kalman  filter  is  only  the  first  iteration  of  the 
single  stage  iterative  estimator. 

Let  xK1  ,K,  ,  denote  the  estimate  of  x..  based 
upon  all  the  measurements  up  to  and  including  .  Upon 

receiving  the  new  measurement  Zjq,  the  optimum  estimate  at 

A 

time  tN,  xN/N/  is  the  xN  minimizing 
J  =  (z-h(xN))TRN_1(zN-h(xN)) 

+  ( -N/N-1  ~  — N^  P N/N-1 ( -N/N-1  -N}  (B,5) 

where  and  PN/N_i  ace  covariances  of  _zN  and  1  '  res_ 

pectively.  Following  similar  derivations  of  [16],  one  obtains 
the  following  iterative  algorithm. 


! 


^k+1  k  lc  —  1  *  k  —1  ~  ~]r 

-N/N  =  -N/N  +  PN/N  HnRn  (-n"-(-N/N) }  +  PN/N- 1 ' — N/N- 1 N/N } 


=  P_1  +  HkR  "1HkT 

N/N  N/N-1  nNN  N 


(B.6) 

(B.7) 


v 

where  is  the  Jacobian  matrix  of  h (  • ) .  Equation  (B.7) 
can  be  re-written  as 


PN/N  PN/N- 1 [ 1  HN  (HNPN/N-1HN  +  ^  )  HNPN/N-1 


(B.8) 


using  the  Matrix  Inversion  Lemma. 

Notice  that  if  we  choose  the  initial  guess  for 

"0 

—N/N  as  —N/N  1  and  stoP  after  the  first  iteration,  we  have 
obtained  the  extended  Kalman  filter  equation. 
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